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ABSTRACT 

A triangle based TVD (total variation diminishing) scheme for the numerical approxima- 
tion of hyperbolic conservation laws in two space dimensions is constructed. The novelty of 
the scheme lies in the nature of the preprocessing of the cell averaged data, which is accom- 
plished via a nearest neighbor linear interpolation followed by a slope limiting procedure. 
Two such limiting procedures are suggested. The resulting method is considerably more 
simple than other triangle based non-oscillatory approximations which, like this scheme, ap- 
proximate the flux up to second order accuracy. Numerical results for linear advection and 
Burgers’ equation are presented. 


1 This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18605 while the second author was in residence at the Institute for Computer Applications 
in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 

a Research supported by Chevron Oil Field Research Company, ONR Grant N00014-86-K-0691, NSF 
Grant DMS 88-11863, DARPA Grant in the ACMP program, and NASA Langley Grant NAG1-270 and 
NASA University Consortium NCA-2372. 


l 



1. Introduction 


In the last ten years there has been considerable effort aimed at constructing 
and analyzing high order accurate, non-oscillatory approximations to hyperbolic 
conservation laws (see e.</., [2], [8]). It is by now well established that the spon- 
taneous development of shock waves and the appearance of steep gradients in the 
solution require higher order schemes to have an adaptive stencil (by adaptive stencil 
we mean an adaptive flux approximation, not an adaptive grid) in order to suppress 
the spurious oscillations that plague conventional finite difference methods. Total 
variation diminishing (TVD) schemes, one such class of second order accurate meth- 
ods that eliminate unphysical oscillations, have been used successfully in a variety 
of applications. Recently, a new class of methods, essentially non-oscillatory (ENO) 
schemes ([3], [4]), which surpass the second order accurate barrier associated with 
TVD schemes, has been developed. An alternative approach for third order schemes 
was developed in [10]. 

Extensions of TVD and ENO schemes to two and three dimensions are typically 
accomplished in a dimension by dimension fashion, via space-operator splitting. 
Therefore, the extension of these higher order schemes to the solution of hyperbolic 
conservation laws on unstructured grids, such as a triangular mesh, is not imme- 
diate. It is our intent in this paper to devise a second order accurate scheme of 
TVD type which is applicable to an unstructured triangular grid. Our scheme is 
based on a finite volume type discretization and is particularly straightforward to 
implement. The scheme relies on a very local adaptive interpolation idea, which 
results in computational efficiency. In the future, we expect to extend the adaptive 
two dimensional interpolation ideas presented here to develop triangle based, higher 
order ENO schemes. 

Several approaches for the solution of hyperbolic conservation laws on trian- 
gular grids already exist. These techniques are, however, in the context of finite 
element methods and have utilized flux corrected transport (FCT) ideas [6] or have 
required the generation of a complex auxiliary grid [9] or are truly finite element 
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methods in space and time and thus are more costly computationally [5], The 
methodology presented here is, in our opinion, simpler and more efficient, primarily 
because a finite volume rather than a finite element approach is used, thus avoiding 
the overhead associated with finite element schemes. 

The TVD, second order accurate methods we shall develop in §2 are technically 
neither total variation diminishing nor strictly second order accurate. We follow the 
convention of calling two dimensional schemes TVD if they are formal extensions 
of one dimensional TVD schemes, as our scheme is. In general, however, the total 
variation may increase [1], though a maximum principle is satisfied. Also, although 
the fluxes are approximated up to second order, the truncation error is technically 
lower because of the adaptive stencil and the variable size of the triangles. Numerical 
experiments presented in §3 indicate orders of accuracy between 1.6 and 1.9 in the 
L\ norm. In §4 we suggest further extensions of the method within the TVD context 
and indicate partial extensions to include diffusive terms. 

2. Construction of the Numerical Schemes 

Our intent in this section is to develop a scheme to solve hyperbolic conserva- 
tion laws on triangular grids in two space dimensions. The method presented is for 
single hyperbolic conservation laws, though hyperbolic systems can be treated anal- 
ogously in a field by field manner. Our method is finite volume based and achieves 
greater than first order accuracy through use of a novel adaptive flux interpolation 
procedure. We first present the general finite volume approach, then introduce our 
general limiting procedure, and then discuss various specific limiters. 

2.1 Finite Volume Discretization 

Consider the hyperbolic conservation law, 

Ut + v- F(u) = g(x,t), 
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u(x,0) = u 0 (x), 


( 2 . 1 ) 


subject to boundary conditions. We wish to solve (2.1) on a triangular grid, a 
portion of which is shown schematically in Fig. 1. Integrating (2.1) over a triangle 
(A abc to be specific) gives, 

^ / udA = - f (V • F)dA, (2.2) 

JAaBC J&ABC 

where A abc represents both the region ABC and its area, and <?(x,t) has been 
taken to be zero for simplicity of exposition only. Applying the divergence theorem 
to the right hand side of (2.2) and defining, 


u = 


( f udA ) (A 

J A arc 


ABC 



i.e., u is the average of u over A abc, gives. 


dt U &ABC 


L 


F • n ab^I + f F • n Acd£, 

L JtAB J^AC 


+ f F • n Bcdi . 
J Irc 


(2.3) 


Note that u is equal to the value of u evaluated at the triangle centroid ( x ^ bc ) 
to within 0{Aabc\ or, analogously, to within 0(^ 2 ), where i is the characteristic 
length of a side of A abc- Here n is the unit outward normal. 


We approximate (2.3) by first using a semi-discrete approach where the ap- 
proximation is 

VABc(t ) ~ UABc(t)\ 

the same is true for all triangles. First order accurate monotone schemes can eas- 
ily be constructed - see e.g., [7], [14]. Let h B c{} 1} \i w 2) be a two-point Lipschitz 
continuous monotone flux, approximating F • n gc, i.e., 

(2.4a) h BC (w,u>) = F • n B c, 

(2.4b) h B c(wi,W 2 ) is a nondecreasing function of w\ and a nonincreasing func- 
tion of U >2 ■ 
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Then our semidiscrete monotone approximation is, 

1 

&ABC 

+ }iab(vabc,vabf) ■ ?ab (2.5) 

+ h,Ac{vABC , VACE ) • f-AC 

where £bc is the length of the side BC, etc. “E” schemes may also be used - see 
[7] for the definition and for examples. 

To obtain higher order accuracy we preprocess our initial data so that in each 
triangle, in particular A ABCi a linear function is obtained whose cell average equals 
vabCi but which is within 0(A) of uabc in regions of smoothness. Here A is the 
maximum area of the four triangles seen in Fig. 1. Moreover, this linear function 
will not introduce new oscillations in our approximation. This (simple) construction 
is the key part of this paper; it will be described at the end of this section. We call 
this linear approximation La( x )- It is generally discontinuous across the boundary 
of each triangle. 

Let ~x.bc be the midpoint of side B0, etc. Let Ba( x BC ) denote the limit of 
La( x ) as x — » Xbc from inside triangle ABC and La( x< bc) denote the limit as 
x — + xbc from outside triangle ABC. Generally, 

|La( x 5 C ) — La(x5c)| = 0(A). 

Our TVD, second order accurate, semi-discrete approximation to (2.3) is 

d 1 ' 

~O^VABc(t) — — ^ ^Bc(La( x bc)> I"a( x Bc)) ■ ^BC 

+ ^ab(La( x !4b)> ^^( x ab)) ‘ ^ab (2-6) 

+ ^ j 4c(La( x ^c)> La(x°ac)) ' &AC • 

By the midpoint formula for integrals, this approximation is weakly second order 
accurate, in the sense that each of the three flux terms above is within 0(A) of 



a?' Msc(<) 


hBc(vABC, VBCD) ' ^BC 
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the line integrals, f F • nd£, along the corresponding interfaces. However due to 
the shifting stencil and varying size and relation of the triangles, the pointwise 
truncation error is generally only 0(A 2 ), i.e., first order. The performance appears 
to be around 1.6- 1.9 order in L\ for smooth flow (see §3). 

2.2 Construction of Linear Function La 


We now describe the construction of La- In each interior triangle, three can- 
didates for La, designated L' A , are generated. The first such candidate is the 
linear interpolate of the three values 

(X-ABC,VABc), (*BCD,VbCd), ( X-ACE,VACE ), 

L\ is the interpolation of 

{xabc,vabc), (xbcd,ubcd), (xabf,vabf), 

and L\ the interpolation of 

C * ABC, V ABC ), ( X-ACE,VACE ), ( X-ABF,V A BF )■ 

These three linear interpolants are sketched in Fig. 2. Here and below we assume 
that the three triangle centroids, ~x.abc ,^bcd ail d xabf are not colinear. At this 
point, three possible L\ exist, and a limited version of La must be selected from 
these. To accomplish this, we first compute the magnitude of the gradient of each 

L a , Le., 



By analogy with limiting procedures in one space dimension ([13]), a valid, though 
very non- compressive limiter, corresponds to the selection of the L\ for which ] V L' A \ 
is the minimum. This choice is analogous to the min limiter in second order ENO 
methods ([3]); no special precautions need be taken at extrema. 
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It is desirable to construct a more compressive limiter than that described 
above, particularly for problems involving linear or contact discontinuities. To ac- 
complish this, we first consider the more compressive slope limiters in one dimension, 
the $ type limiters described by Sweby [13] (his equation 3.17) of which superbee 
is the most compressive, corresponding to $ = 2. These limiters allow the use of 
piecewise linear approximations to the solution for which the slope is not the min- 
imum, subject to the restriction that no overshoot (or undershoot) occurs at the 
cell boundaries. 


The next limiter we describe is a multidimension analog of the one dimensional 
$ limiters. The approach here is to select the L A for which |VZ^| is maximized, 
subject to the restriction that no overshoot or undershoot occurs at any of the three 
triangle boundaries. The procedure is as follows: 

(i) Select the L A for which |VZ^| is the maximum. 

(ii) Check for overshoot or undershoot at ~x.ab,*ac and xgc 1 - For L l A to represent 
a valid Za, it suffices to verify that, for A abc, 

La{*ac) is between vabc and vace. 

La (xab) is between vabc and vabf and 

Ta(xbc) is between vabc and vbcd- 

If these three requirements are satisfied, L' A is the appropriate La 

(iii) If the L\ above results in overshoot or undershoot at any one of the three 
midpoints, select the L A for which |VZ^J is the second largest and repeat the 
test in (ii). If this L l A does not satisfy the test in (ii), select the L A for which 
|VZ^| is the minimum and again proceed through the test in (ii). 

(iv) If all L a fail (ii), revert to a piecewise constant approximation for A abc', i.e., 
La = vabc- 

Given Za, the right hand side of (2.6) can be evaluated and vabc (f) integrated 
in time. This time integration is accomplished via a second order TVD Runge-Kutta 
procedure [11]. 
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3. Numerical Verification of Higher Order Scheme 


In this section we present results for the convergence of the general method 
described in §2, as well as solution contours and profiles demonstrating the accuracy 
of the method. In all cases, the solution region is a square domain discretized via 
right triangular ‘volumes’ (referred to as elements), as shown in Fig. 3. Periodic 
boundary conditions are imposed in both the x— and the y— directions; the initial 
condition is similarly x— and y— periodic. In all cases the more compressive limiter 
described in §2.2 is used for the TVD scheme. 

S.l Rate of Convergence 

To assess rate of convergence, the scheme is applied to the solution of the linear 


conservation law 


Ut -p V • (au) = 0, 

(3.1) 

subject to the initial condition 


u 0 (x,y) = sin(27rx)sin(27ry). 

(3.2) 

Our base monotone scheme uses the EO flux [7]: 


h(wi,w 2 ) = f+(w i) + f-(w 2 ). 

(3.3) 

For linear equations with constant a = (a x ,a y )^ 


f+(u ) = [max((a • n),0)]u, 

(3.4a) 

/_(u) = [min((a • n), 0)]«. 

(3.46) 


For Burgers’ equation (considered below), where /i = /2 = (l/2)t/ 2 , 

/+(«) = rciax(( n * ^ — )u 2 , 0), (3.5a) 

/-(«) = min(( — + - -- )u 2 ,0), (3.56) 
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where n x and n y represent the components of n. 


A contour plot of the initial condition (3.2) is shown in Fig. 4. Four extrema 
are evident. The rate of convergence of the method was determined for both the 
case a x = a y = 1 and a x = 1, a y = 0. Further, convergence was assessed both on 
an element by element basis and after applying a local averaging procedure. It is 
expected that local averaging procedures would enhance the rate of convergence, as 
the scheme is expected to be second order in only the weak sense; i.e ., after integrat- 
ing locally in space and time (a type of local averaging). The averaging performed 
in this study is, however, only spatial; no temporal averaging is attempted. This 
is because spatial-temporal averages are rather cumbersome to perform in practice, 
and the spatial averaging alone reveals the expected trend. Computations were 
performed for grids ranging in discretization from 200 elements (£ = 0.1, where £ 
is the spacing between adjacent nodes or, analogously, £ = (2A) 1 / 2 , with A the 
area of any element) to 12800 elements {£ = 0.0125). In all cases the CFL number, 
A(= At /£), was set to 0.1. 

Displayed in Fig. 5 is a log-log plot of L\ error versus £. In this case, a x = a y = 
1. Results are shown for both a first order scheme and the higher order scheme, 
with error computed on an element by element basis. Least squares linear fits give 
the order of convergence for the two methods; for the first order method we obtain 
0.93 and for the higher order method 1.77. Figure 6 displays an analogous plot 
after applying a local averaging procedure. Specifically, this averaging procedure 
entails averaging the computed value of u over square regions comprised of two 
adjacent elements and computing the error in terms of the difference between this 
average and the exact solution of Eq. (3.1) evaluated at the square midpoint. For 
the grid displayed in Fig. 3, 100 such square regions exist. Again, averaging is 
only applied spatially; no temporal averaging is performed. Assessing error in this 
manner results in least squares linear fits of slope 0.94 for the first order method and 
1.81 for the higher order method. As expected, local averaging enhances the rate 
of convergence though, in this case, the improvement is minimal. In other cases, 
however, the improvement is more substantial. For example, using a x — 1, a y = 0 
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in Eq. (3.1) yields the following convergence results. For the first order method, 
convergence is 0(0.99) with no local averaging and 0(1.06) with averaging. For the 
higher order scheme, the convergence rates are 0(1.60) and 0(1.73), respectively. 

Results for L 2 and error display slower rates of convergence. For the case 
a x = a y = 1, L 2 error is ~ 0(£ 164 ) with local averaging and 0(£ 162 ) on an element 
by element basis with Loo error ~ 0(£°' 94 ) with local averaging and 0(£°' 91 ) with no 
averaging. The expected result for L 2 error is ~ ()(£ l ' 5 ) and for Loo error ~ 0(7), as 
in one dimensional TVD methods. Although the discrepancies between the expected 
and numerical results are relatively slight, it is not clear why the L 2 error converges 
faster than expected while the L 0 0 error converges slower than expected. 

Shown in Table 1 is a compilation of the rates of convergence of L \ , L 2 and 
error. Results for both a first order scheme and our more compressive TVD scheme 
are displayed. In all cases the initial condition is as in (3.2). Error is computed 
over the entire domain in two ways: (1) element by element and (2) by combining 
two adjacent elements into squares. In all cases i ranges from 0.0125 to 0.1. 

Slightly improved rates of convergence in L\ are obtained when the initial 
condition contains no extrema. This is demonstrated in Table 2, where results 
for Li error for the initial condition Uo(x,y ) = sin(7ra:/2)sin(7ry/2) are displayed. 
Here, to eliminate the effects of the discontinuity in u at the boundary (recall that 
periodic boundary conditions are imposed), error is computed only over the region 
0.6 < x, y < 0.8 at an early time, t = 0.05. In one case, local averaging has a more 
dramatic effect, improving the L\ accuracy of the TVD scheme from 1.22 to 1.80. 

Based on the numerical results presented above and the analysis presented in 
§2, we feel that the method can be considered to be second order accurate in L\ in 
the weak sense. Though our convergence results always indicate convergence slower 
than quadratic, this is, in our opinion, due to the fact that these results are not 
strictly measuring weak convergence. If such convergence could be unambiguously 
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measured, it is our contention that the method would indeed display second order 
convergence. 

3.2 Examples of Numerical Accuracy 

We now present some detailed numerical results for our second order scheme 
and compare these with the results of a first order method. The first results are for 
the solution of Eqs. (3.1) and (3.2) with a x = ay = 1. Figure 7 displays the solution 
contour results for the first order scheme with 800 elements (£ — 0.05) and A = 0.1 
(the same CFL number is used in all computations) at t = 1. The exact solution 
is a reproduction of the initial condition, shown in Fig. 4. The first order method 
is clearly very diffusive; the maximum value of u is here only 0.25, in contrast to 
the maximum in the initial condition of 1. Results for the second order scheme 
at t = 1 are shown in Fig. 8. Though some distortion of the initial condition is 
apparent, the solution is considerably improved over the first order solution; the 
maximum value of u is now 0.76. Shown in Fig. 9 are the t = 1 results for the 
first order scheme using 3200 elements (£ = 0.025). Substantial numerical diffusion 
is still evident; the maximum value of u is only 0.49. The solution contour using 
the second order method is displayed in Fig. 10. The t = 1 solution in this case 
closely resembles the initial condition, with a maximum value of u of 0.88. Figures 
11 and 12 show solution profiles taken along the line y — x (the velocity direction) 
at t = 0, 0.25, 0.5, 0.75 and 1 for both the first and second order methods. In 

both cases, 3200 elements were used. The second order results are quite sharp at 
all times, while the first order results show a continual degradation with increasing 
time. 

Solution profiles for computations using the second order method with a non- 
linear flux function, / = (l/2)u 2 in Eq. (2.1) (be., the inviscid Burgers’ equation), 
with the initial condition uo(x,y) = sin(27rx), are shown in Fig. 13. Though this 
is an essentially one dimensional problem, no overshoot or unphysical oscillations 
appear in the solution. The solution is sharper with the second order method than 
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with a first order scheme (first order results are not shown), though the improvement 
is of course less dramatic than in the linear examples presented above. 


4. Possible Extensions 


Other limiting procedures are quite feasible and should be tested. Our com- 
pressive limiter is not a direct analogue of superbee, since superbee (and many other 
limiters [13]) occasionally allows values other than zero or any of the slopes being 
compared to be the final choice of slope (or gradient in our two dimensional case). 


A more significant issue is the treatment of diffusive terms. In this case, the 
governing equation is of the form 


u t 


&ABC 


F(u) = e(u xx + u yy ), e > 0. 

(4.1) 

now involves the additional term, 


fi de+ 1 %-dt+f 

dn J tAC dn J tBC dn 

(4.2) 


on the right side of (2.3). Up to first order accuracy, we compute each of the three 
terms in (4.2) as follows. The limiting procedure has already given us a gradient 
within the triangle ABC as well as for each of the three neighbors. Therefore, the 
integral along side AB in (4.2) can be computed approximately as 


[(VLabc + V Labf ) • 2 ~ ■ 


(4.3) 


The integrals along the other sides are approximated analogously. This is generally 
a first order accurate method (second order accuracy occurs in special cases, t.g., if 
all the triangles are equilateral). However, since e is relatively small here (otherwise 
transport is diffusion dominated and the sophisticated treatment of convection is 
unnecessary), we believe this to be an adequate treatment of these terms. 


Finally, we mention that work is underway to approximate (2.1) using a higher 
order accurate ENO triangle based method. See [11], [12] for successful Cartesian 
coordinate approaches. 
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Table 1 . Computed accuracy of TVD scheme for the linear case. Initial condition 
u o( x >y) — sin(27ra:) sin(27ry). Error computed over the entire domain. 
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Table 2. Computed accuracy (L\) of TVD scheme for the linear case. Initial 
condition uo(x,y) = sin(7rx/2) sin(7ry/2) contains no extrema. Error computed 
over 0.6 < x, y < 0.8 at t = 0.05. 
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Figure 3 

Triangular grid used for the numerical calculations. 
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Figure 4 

Contour plot of the initial condition (3.2). Contours correspond 
to u = 0, ±0.15, ±0.3, ±0.45, ±0.6, ±0.75, ±0.9. 
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Figure 5 

L i error on a per element basis for the case a x = a y = 1 1 
first order (O) and second order (x) schemes. Lines are 
least square fits with slopes as indicated. 


Li error after local averaging for the case a x = a y = 1 foi 
first order (O) and second order (>$ schemes. Lines are 
least square fits with slopes as indicated. 
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Figure 1 1 

Solution profiles along the line y=x for the first order scheme 
(3200 elements). 
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Figure 13 

Solution profiles for Burgers' equation using second order 
scheme (800 elements). 
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